#define _IMPCAPT_
!#define RRTALLY
module util_STT
    implicit none
contains

function var(buf, st, ed)
    real(8), intent(in), dimension(:) :: buf
    integer, intent(in) :: st, ed
    integer :: i,n
    real(8) :: sqsum, sum, var
    n=ed-st+1
    sqsum=0
    sum=0
    do i=st, ed
        sqsum = sqsum + buf(i)**2
        sum = sum+buf(i)
    end do
    var = (sqsum-sum**2/n)/(n-1)
end function

function stdev(buf, st, ed)
    real(8), intent(in), dimension(:) :: buf
    integer, intent(in) :: st, ed
    integer :: i,n
    real(8) :: var, sqsum, sum, stdev
    n=ed-st+1
    sqsum=0
    sum=0
    do i=st, ed
        sqsum = sqsum + buf(i)**2
        sum = sum+buf(i)
    end do
    var = (sqsum-sum**2/n)/(n-1)
    stdev=sqrt(var)
end function

function stdevp(buf, st, ed)
    real(8), intent(in), dimension(:) :: buf
    integer, intent(in) :: st, ed
    integer :: i,n
    real(8) :: var, sqsum, sum, stdevp
    n=ed-st+1
    sqsum=0
    sum=0
    do i=st, ed
        sqsum = sqsum + buf(i)**2
        sum = sum+buf(i)
    end do
    var = (sqsum-sum**2/n)/(n)
    stdevp=sqrt(var)
end function

function covarlag(buf, st, ed, lag)
    real(8), intent(in), dimension(:) :: buf
    integer, intent(in) :: st, ed, lag
    real(8) :: xt, xs, psum, covarlag
    integer :: i, n
    n = ed-st-lag+1
    xt = sum(buf(st : ed-lag))/n
    xs = sum(buf(st+lag : ed))/n
    psum = 0.
    do i=st, ed-lag
        psum = psum+buf(i)*buf(i+lag)
    end do
    psum = psum/n
    covarlag = psum-xt*xs
end function

end module

module util_mc
    use define
    use RNG
    implicit none
    integer, private, parameter :: fphinm = 101, fsptr = 71
!    integer, private, parameter :: maxring=2
!    real(8), private, dimension(maxring) :: rfuel, rfsq
    integer, private, parameter :: rin=1 , rout=2
    integer, private, parameter :: PSINTO=1, PSSTAY=2, PSOUT=3
    integer, private, parameter :: INDIR=1, OUTDIR=2
    integer, parameter :: TTPASS =0
    integer, parameter :: TTREFL =1
    real(8), pointer, dimension(:,:,:) :: psicm
!    integer, save :: testiter = 0
contains

subroutine init_util_mc
    use allocs
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    integer :: ig, ic
!    rfuel(1) = hpx * 0.8 / 2    
!    rfuel(2) = hpx * 0.6 / 2    
!    rfsq = rfuel**2
    open (fphinm, file = "phinm.txt", status = 'unknown')
    open (fsptr,  file = "sptr.txt", status = 'unknown')
    call dmalloc(psicm, ncmx, ncmy, 1)
end subroutine

subroutine migrate(neutron, k, ke, nnt, bseed, mi &
!$    ,tid, tlpsi, tljc, tlk, tlntrn, tlnq &

    )
    use asmcal
    implicit none                                                              
    include "arrays.FI"    
!    include "mtbuf.FI"                                                    
    include "pspec.FI"                                                         
    include "mcbase.FI"  
    include "xsec.FI"    
    include "FMC.FI"         
    include "GroupConstant.FI"                              
    real(8), parameter :: mu_min = 0.01_8                                                                               
    type(particle), intent(in) :: neutron                                      
    real(8), intent(in) :: k, ke                                             
    integer, intent(out) :: nnt  
    integer(8), intent(in) :: bseed
    integer, intent(in) :: mi
!    integer(1) :: l, m ! , lm(2)
    integer, dimension(2,4) :: LTRFLT
    data LTRFLT /-1,1, -1,1, 1,-1, 1,-1/
                                             
    logical, parameter :: ABSR_ = .true.                                        
    real(8) :: dist, dist_surf, alpha, trkl, test                          
    integer :: i, gtilde, surf, btype, nn, type_col, isurf, ttt, inc1, inc2 
    integer :: nnw, ifmc
    logical escape          
    real(8), parameter :: p=0.5, s= 0.2
    real(8) :: sigts, mud
    logical :: fly
    real(8) :: phixx, phixy, phiyy
    real(8) :: phs
    ! equation for neutron trk line
    real(8) :: a, b, x0, x1, nwt
    integer(8) :: lseed

!    ! index for coarse mesh surface (Larsen's CMFD)
!    integer :: si, sj, axis
!    real(8) :: wphi2, wtr
!
    ! module -> subroutine
    real(8) :: xy(2), x, y, dir(2), mu, pr, w, mst, at
    integer :: lm(2), g, tcomp, l,m, lc, mc , gc
    integer :: nloop

    integer :: inout, coltype, ri, tpin
    real(8) :: dts, rx, ry, ox, oy, rsq
    integer :: psts ! position status
    integer :: rd
    logical :: terminate
    type (strcpin), pointer :: ppin
    real(8), pointer, dimension(:,:,:) :: rtot, rnuf, rabsr, phiv
    real(8), pointer, dimension(:,:,:,:) :: rsct
    real(8), pointer, dimension(:,:,:) :: pau, psuu
    real(8), pointer, dimension(:,:,:) :: prcwx, prcwy, pphi2x, pphi2y, ptrx, ptry
    logical :: lleak
!$  integer, intent(in) :: tid
!!$  integer(8), intent(in) :: lseed
!!$  real(8), dimension(0:nmxr,1:ng,1:nfmx,1:nfmy)      :: tltrk
!$  real(8), dimension(nfmx, nfmy)                      :: tlpsi
!$  real(8), dimension(1:4, 1:ngcmfd, 0:ncmx+1, 0:ncmy+1)    :: tljc
!$  real(8)                          :: tlk
!$  type(particle),dimension(:,:)    :: tlntrn
!$  integer, dimension(:)            :: tlnq

    real(8), pointer, dimension(:,:,:) :: asmphiv, asmtot, asmnuf, asmabs
    real(8), pointer, dimension(:,:,:,:) :: amssct
    integer :: gf
    

    equivalence (l, lm(1))
    equivalence (m, lm(2))
    equivalence (x, xy(1))
    equivalence (y, xy(2))
    equivalence (ox, dir(1))
    equivalence (oy, dir(2))
    
    lleak = .false.
    

    if (tallytype .eq. MD_MCCMFD)then
        if (hmzt .eq. HMZMC) then
            rtot=>GCTtot(:,:,:,tid)
            rnuf=>GCTnuf(:,:,:,tid)
            rabsr=>GCTabsr(:,:,:,tid)
            phiv=>GCTphiv(:,:,:,tid)
            rsct=>GCTsct(:,:,:,:,tid)
            pau=>tphiau(:,:,:,tid)
            psuu=>tphisuu(:,:,:,tid)
            
            ! GCMFD
            prcwx=>ttrcwx(:,:,:,tid)
            prcwy=>ttrcwy(:,:,:,tid)
            pphi2x=>tphi2x(:,:,:,tid)
            pphi2y=>tphi2y(:,:,:,tid)
            ptrx=>tphi2trx(:,:,:,tid)
            ptry=>tphi2try(:,:,:,tid)
        endif
    elseif (tallytype .eq. MD_ASMCAL)then 
        call sync_asmarray(tid, asmphiv, asmtot, asmnuf, asmabs, amssct)
    else
        
    endif
!    equivalence (l, lm(1))
!    equivalence (m, lm(2)) 


!   Initialization of fission source before migration ====================================
    xy = neutron%xy
    w = neutron%w                                                              
    lm = neutron%lm 

    rx = x -(hacx(l) - hfmx/2)
    ry = y -(hacy(m) - hfmy/2)
    rsq = sqrt(rx**2+ry**2)
    tpin=cplst(l,m)
    ppin=>typepin(tpin)

    ! check the pstatus of source neutron
    do ri=ppin%nr, 1, -1
        inout = inout_(ri)
        if (inout .eq. rin) exit
    enddo

    tcomp=ppin%rc(ri)

    ! skip stride
    lseed=strideN(bseed, mi)
    g=nn_group()
    mu = 1-2*GetRNS(lseed)
    mst = sqrt(1._8-mu**2)
    at = 2*pi*GetRNS(lseed)
    dir(1) = mst*cos(at)
    dir(2) = mst*sin(at)

    lc=f2cx(l)
    mc=f2cy(m)
    gc=gcmc(g)

    nnt = 0                                                                    
    escape = false 

#ifdef ASMPRL            
    tnfs(lc,mc,tid)=tnfs(lc,mc,tid)+1
#endif
!=========================================================================================


    dist = -1._8/Xtot(tcomp, g)  * log(GetRNS(lseed))   

    terminate = .false.
    rd = INDIR
    
    do while (terminate .eq. .false.)
        if (rd .eq. INDIR) then
            dts = dtcl_(ri+1, rout)
            if (dts .eq. 0._8)then 
                rd = OUTDIR
                continue
            else
                ! inside
                if (dist >dts)then
                    ! move to inner ring
                    dist = dist -dts
                    call move(dts)
                    ri = ri+1
                    call check_material
                    cycle
                else
                    ! collision
                    call move(dist)
                    coltype = collision()
                    if (coltype .eq. ABSR_) terminate = .true.
                    cycle
                endif
            endif
        endif
        if (rd .eq. OUTDIR) then 
            ! outside 
            if (ri.eq.0)then; dts=dist2surf(surf)
            else;dts = dtcl_(ri, rin)
            endif
            if (dist>dts)then
                ! move to outer ring
                dist = dist -dts
                call move(dts)
                if (ri.eq.0) then; 
                    terminate=movenext(surf);
                    rd=INDIR;
                else; ri = ri-1;
                    call check_material
                endif
                cycle
            else
                ! collision
                call move(dist)
                coltype = collision()
                if (coltype .eq. ABSR_) then
                    terminate = .true.
!                else
!                    call check_material
                endif
                rd=INDIR
                cycle
            endif
        endif
    enddo
contains

subroutine move(dst)
    include "GroupConstant.FI"
    real(8), intent(in) :: dst
    real(8) :: trkl
    integer :: gp, gf2
    xy = xy + dst*dir
    rx = x -(hacx(l) - hfmx/2)
    ry = y -(hacy(m) - hfmy/2)
    rsq = sqrt(rx**2+ry**2)
    
    trkl = dst*w
!    if (trkl<0) then
!        print *, "181818"
!    endif
#ifdef _OPENMP
!    tltrk(ri,g,l,m) = tltrk(ri,g,l,m) + trkl
    if (tallytype .eq. MD_MCCMFD) then
        tlpsi(l,m) = tlpsi(l,m) + Xnuf(tcomp,g)*trkl
#ifndef _NOCMFD    
        if (hmzt .eq. HMZMC) then
            phiv(gc,lc,mc)   =phiv(gc,lc,mc)   +trkl
            rtot(gc,lc,mc)   =rtot(gc,lc,mc)   +Xtot(tcomp,g)*trkl
            rnuf(gc,lc,mc)   =rnuf(gc,lc,mc)   +Xnuf(tcomp,g)*trkl
            rabsr(gc,lc,mc)  =rabsr(gc,lc,mc)  +Xabs(tcomp,g)*trkl
            if(lmgcmfd)then
                do gp=1, ng
                rsct(gp,gc,lc,mc) =rsct(gp,gc,lc,mc) +xsctm(gp,g,tcomp)*trkl
                enddo
            else
                rsct(1,gc,lc,mc) =rsct(1,gc,lc,mc) +GCsctm2g(1,g,tcomp)*trkl
                rsct(2,gc,lc,mc) =rsct(2,gc,lc,mc) +GCsctm2g(2,g,tcomp)*trkl
            endif
            
            call track_volume_tallies_GCMFD(trkl)
        endif
#endif    
#else
        trk_est(l,m,g) = trk_est(l,m, g) + trkl  
#endif
    elseif(tallytype .eq. MD_ASMCAL) then
        gf=gtog(g)
        asmphiv(gf,l,m)=asmphiv(gf,l,m)  +trkl
        asmtot(gf,l,m) =asmtot(gf,l,m)   +Xtot(tcomp,g)*trkl
        asmnuf(gf,l,m) =asmnuf(gf,l,m)   +Xnuf(tcomp,g)*trkl
        asmabs(gf,l,m) =asmabs(gf,l,m)   +Xabs(tcomp,g)*trkl
        do gp=1, ng
            gf2=gtog(gp)
            amssct(gf2,gf,l,m) =amssct(gf2,gf,l,m) +xsctm(gp,g,tcomp)*trkl
        enddo
    else
    endif
endsubroutine


subroutine CENTREFL(surf)
    integer, intent(inout) :: surf
    integer :: REFLSURF1(4)
    integer :: REFLSURF2(4)
    data REFLSURF1 / DE, DE, DN, DS/
    data REFLSURF2 / DW, DE, DS, DN/
    select case (reflpin(l,m))
    case(1)
        if (rx <0) then
            surf=REFLSURF1(surf)
            rx=rx*-1.0_8
            ox=ox*-1.0_8
            x=x+rx*2
        endif
    case(2)
        if (ry >0) then
            surf=REFLSURF2(surf)
            ry=ry*-1.0_8
            oy=oy*-1.0_8
            y=y+ry*2
        endif
    case(3)
        if (rx <0) then
            surf=REFLSURF1(surf)
            rx=rx*-1.0_8
            ox=ox*-1.0_8
            x=x+rx*2
        endif
        if (ry >0) then
            surf=REFLSURF2(surf)
            ry=ry*-1.0_8
            oy=oy*-1.0_8
            y=y+ry*2
        endif
    end select
end subroutine

function movenext(surf) result (term)
    include "pspec.FI"
    integer, intent(inout) :: surf
    logical :: term
    integer :: tlc, tmc
    term = FALSE

!#ifdef _OPENMP            
!    jfomp(surf,l,m,g)=jfomp(surf,l,m,g)+w    
!#else
!    jf(surf,l,m, g)=jf(surf,l,m, g)+w    
!#endif


    ! treatment for CENT reflection
!    if ( refltype .eq. CENT) then
    if (refltype_pin .eq. CENT)then
        if (reflpin(l,m)>0) call CENTREFL(surf)
    endif
!    endif
    
    ! boundary check
    if (chk_bndry(surf)) then ! boundary
        if (albedo(surf)<0._8)  then ! PERIODIC================================================
            call move_periodic(surf)
            tpin=cplst(l,m)
            ppin=>typepin(tpin)
            call check_material                           
            if (tallytype .eq. MD_MCCMFD)then
#ifndef _NOCMFD ! ???
                if( (lc.ne.f2cx(l)) .or. (mc.ne.f2cy(m)) ) then ! if the coarse mesh index was changed
                    tljc(surf,gc,lc,mc)=tljc(surf,gc,lc,mc)+w
#ifdef ASMPRL   
                    if ( .not.(lleak)) then         
                        tnout(lc, mc,tid)=tnout(lc, mc,tid)+1
                        lleak = .true.
                    endif
#endif
                    ! larsen's CMFD
                    call track_surface_tallies_GCMFD(TTPASS)
                endif
#endif
            endif
        elseif (albedo(surf)<1e-10_8) then  ! REFLECTIVE================================================

            if (tallytype .eq. MD_MCCMFD) call track_surface_tallies_GCMFD(TTREFL)

            dir = dir * LTRFLT(:,surf)                                             
        else ! ESCAPE ================================================
            ! escape
            term = TRUE                                              
#ifndef NOCMFD ! NO CHANGE IN COARSE MESH INDEX   
            tljc(surf,gc,lc,mc)=tljc(surf,gc,lc,mc)+w
#ifdef ASMPRL            
!            tnout(lc, mc,tid)=tnout(lc, mc,tid)+1
            if ( .not.(lleak)) then         
                tnout(lc, mc,tid)=tnout(lc, mc,tid)+1
                lleak = .true.
            endif
#endif
            
            ! larsen's CMFD
            if (tallytype .eq. MD_MCCMFD) call track_surface_tallies_GCMFD(TTPASS)
            
#endif
        end if                                                         
    else ! inner mesh ================================================
        lm = lm + LTNEIGH(:,surf)  
        tpin=cplst(l,m)
        if (tpin .eq. 0) then 
            term = TRUE
        else
            ppin=>typepin(tpin)
            call check_material                           
        endif
        if (tallytype .eq. MD_MCCMFD)then
#ifndef _NOCMFD 
        if( (lc.ne.f2cx(l)) .or. (mc.ne.f2cy(m)) ) then
            tljc(surf,gc,lc,mc)=tljc(surf,gc,lc,mc)+w
#ifdef ASMPRL  
            ! leaving          
!            tnout(lc, mc,tid)=tnout(lc, mc,tid)+1
            if ( .not.(lleak)) then         
                tnout(lc, mc,tid)=tnout(lc, mc,tid)+1
                lleak = .true.
            endif
#endif
            ! larsen's CMFD 
            call track_surface_tallies_GCMFD(TTPASS)

            lc=f2cx(l)
            mc=f2cy(m)
#ifdef ASMPRL 
            ! incoming           
            if (term .eq. .false.) tnin(lc, mc,tid)=tnin(lc, mc,tid)+1
#endif
        endif
#endif
        endif
    endif                            
end function

subroutine track_volume_tallies_GCMFD(trkl)
    real(8), intent(in) :: trkl
    prcwx(gc,lc,mc)=prcwx(gc,lc,mc)+Xtot(tcomp,g)*trkl*ox
    prcwy(gc,lc,mc)=prcwy(gc,lc,mc)+Xtot(tcomp,g)*trkl*oy
end subroutine


subroutine track_surface_tallies_GCMFD (tallytype)
    integer :: tallytype
    integer :: si, sj, axis
    real(8) :: wphi2, wtr
    
    si=lc+1;sj=mc+1
    if (surf>2) then
        sj=sj+(LTSIDXY(surf)-1)/2
    else
        si=si+(LTSIDXX(surf)-1)/2
    endif
!    wphi2=w*(LTPHI2X(surf)*ox+LTPHI2Y(surf)*oy)
!    wtr =w*(LTPHI2X(surf)*oy+LTPHI2Y(surf)*ox)
    
    if (tallytype .eq. TTREFL) then
        wphi2=wphi2*2
        wtr=0._8
    endif
    
    if(LTSURFDIR(surf) .eq. 1) then ! x direction
!        wphi2=w*LTPHI2X(surf)*ox
!        wtr  =w*LTPHI2X(surf)*oy
        wphi2=w/abs(ox)*ox*ox
        wtr  =w/abs(ox)*ox*oy
        pphi2x(gc,si,sj)=pphi2x(gc,si,sj)+wphi2
        ptry(gc,si,sj)=ptry(gc,si,sj)+wtr
    else                            ! y direction
!        wphi2=w*LTPHI2Y(surf)*oy
!        wtr  =w*LTPHI2Y(surf)*ox
        wphi2=w/abs(oy)*oy*oy
        wtr  =w/abs(oy)*oy*ox
        pphi2y(gc,si,sj)=pphi2y(gc,si,sj)+wphi2
        ptrx(gc,si,sj)=ptrx(gc,si,sj)+wtr
    endif
endsubroutine


function dtcl_(ri, inout) 
    integer, intent(in) :: ri, inout
    real(8) :: sgn(2)
    data sgn / +1.0_8, -1.0_8/
    real(8) :: a, b, c, d, t
    real(8) :: dtcl_, dtcmp
    dtcl_ = 0._8
    if (ri>ppin%nr) return
    rx = x -(hacx(l) - hfmx/2)
    ry = y -(hacy(m) - hfmy/2)
    rsq = sqrt(rx**2+ry**2)
    a = ox*ox+oy*oy
    b = (ox*rx+oy*ry)
    c = rx*rx+ry*ry-ppin%rr(ri)**2
    d = b*b-a*c
    if (d>0._8) then
        t=(sgn(inout)*sqrt(d)-b)/a
        if (inout .eq. rin .or. b < 0) dtcl_ = sqrt((ox*t)**2+(oy*t)**2+(mu*t)**2)
    endif
endfunction

function inout_(ri)
    integer :: inout_
    integer, intent(in) :: ri
    inout_ = rout
    if (ri>ppin%nr) return
    if (rx*rx+ry*ry <ppin%rr(ri)**2) inout_ = rin 
endfunction

function collision() result(labsr)
    logical :: labsr
    integer :: gct
    labsr=.false.
    
! tally for 2G CMFD ================================
#ifdef RRTALLY    
    rrtot(gc,lc,mc)=rrtot(gc,lc,mc)+w
    hmphit(gc,lc,mc)=hmphit(gc,lc,mc)+Xtotr(tcomp,g)*w
    rrnuf(gc,lc,mc)=rrnuf(gc,lc,mc)+Xnfot(tcomp,g)*w
#endif
! ==================================================
    
    call sample_Nneutron(k, ke, nn, nnw, nwt)  
    nnt = nnw+nn  
    if (nnt>0) then   
        do i=1, nn                                                         
                call addntrn(OMPNXT)                                
        end do     
        do i=1, nnw                                                        
            call addntrn(OMPCUR)                                
        end do    
    endif
#ifdef _IMPCAPT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    w = w*Xsct(tcomp, g)/Xtot(tcomp, g)
    ! Russian roullete
    if (w < 0.25_8) then
#ifdef _OPENMP
        if(GetRNS(lseed)<w/0.5_8) then
#else
        if(GetRN()<w/0.5_8) then
#endif
            w = 0.5_8
        else
            w = 0.
            labsr = .true.
            return
        endif
    endif
    gtilde = sample_scattering_group() 
    gct=gcmc(gtilde)
    g = gtilde
    gc=gct
    call sample_direction
#ifdef _OPENMP        
    dist = -1._8/Xtot(tcomp, g)  * log(GetRNS(lseed))   
#else
    dist = -1._8/Xtot(tcomp, g)  * log(GetRN())   
#endif  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#else            
    type_col = sample_collision()                                
    if (type_col.eq.COL_AB) then 
    
#ifdef ASMPRL 
            ! incoming           
            tnabsr(lc, mc,tid)=tnabsr(lc, mc,tid)+1
#endif
       
        labsr = .true.
!! tally for 2G CMFD ================================
!#ifdef RRTALLY    
!        rrabs(gc,lc,mc)=rrabs(gc,lc,mc)+w
!#endif
! ==================================================
    else                                                               
        gtilde = sample_scattering_group()                             
        gct=gcmc(gtilde)
!! tally for 2G CMFD ================================
!#ifdef RRTALLY    
!        rrsct(gct, gc,lc,mc)=rrsct(gct, gc, lc,mc)+w
!#endif
!! ==================================================
        g = gtilde
        gc=gct
        call sample_direction
#ifdef _OPENMP        
        dist = -1._8/Xtot(tcomp, g)  * log(GetRNS(lseed))   
#else
        dist = -1._8/Xtot(tcomp, g)  * log(GetRN())   
#endif  
    end if                                                             
#endif
    
endfunction

subroutine check_material
    integer :: comptypen
    ! tpin, ri ! from outer subroutine
    comptypen=typepin(tpin)%rc(ri)
    if (tcomp .ne. comptypen) then                                
        tcomp = comptypen                                       
#ifdef _OPENMP
        dist = -1._8/Xtot(tcomp, g) * log(GetRNS(lseed))                                 
#else
        dist = -1._8/Xtot(tcomp, g) * log(GetRN())                                 
#endif
    end if                                                         
endsubroutine                                                                  

subroutine move_periodic(surf)
    include "pspec.FI"
    integer, intent(in) :: surf
    select case(surf)
    case(DN);m=1;y=0._8
    case(DS);m=nfmy;y=hacy(m)
    case(DW);l=nfmx;x=hacx(l)
    case(DE);l=1;x=0._8
    end select
end subroutine
  
subroutine sample_direction
    implicit none
    real(8) :: a, pr
#ifdef _OPENMP
    mu = 1._8- 2._8* GetRNS(lseed)
    a  = 2._8*pi*GetRNS(lseed)
#else
    mu = 1._8- 2._8* GetRN()
    a  = 2._8*pi*GetRN()
#endif
    pr = sqrt(1._8-mu**2)
    dir(1) = pr*cos(a)
    dir(2) = pr*sin(a)
end subroutine

function chk_bndry(surf)
    include "pspec.FI"
    integer :: surf, lmn(2)
    logical :: chk_bndry
    lmn=lm+LTNEIGH(:,surf)
    chk_bndry = FALSE
    select case (surf)
        case (DN)
            if (lmn(2) > ibndry(surf)) chk_bndry = TRUE
        case (DS)
            if (lmn(2) < ibndry(surf)) chk_bndry = TRUE
        case (DW)
            if (lmn(1) < ibndry(surf)) chk_bndry = TRUE
        case (DE)
            if (lmn(1) > ibndry(surf)) chk_bndry = TRUE
    end select
end function 

function dist2surf(surf) result(ret)
    use define
    implicit none
    include "pspec.FI"
    integer,intent(inout) :: surf
    real(8) :: ret
    integer :: mi, mj, inear
    real(8) :: lxy(2)
    integer :: LTCHK(2,4)
    data LTCHK /DW,DS, DW,DN, DE,DS, DE,DN/
    integer :: quad
    mi=-1; mj=-1;
    if (dir(1)>0)mi=0
    if (dir(2)>0)mj=0
    lxy(1) = (hacx(l+mi)-x)/dir(1)
    lxy(2) = (hacy(m+mj)-y)/dir(2)
    inear = 1
    if (lxy(1)>lxy(2)) inear = 2
    surf = LTCHK(inear, 2*mi+mj+4)
    ret = lxy(inear)
end function

function sample_collision()
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    integer :: sample_collision
    real(8) xi, tot, sct
#ifdef _OPENMP
    xi  = GetRNS(lseed)
#else
    xi  = GetRN()
#endif    
    tot = Xtot(tcomp, g)
    sct = Xsct(tcomp, g)
    if (xi>sct/tot) then
        sample_collision = COL_AB
    else 
        sample_collision = COL_SC
    end if 
end function 

integer function sample_scattering_group
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    integer gp
    real(8) sum, xi
#ifdef _OPENMP
    xi  = GetRNS(lseed)
#else
    xi  = GetRN()
#endif    
!    sum = 0.
!    do gp=1, ng-1
!        sum = sum + XsctM(gp, g, tcomp)/Xsct(tcomp, g)
!        if (xi<sum) then
!            exit
!        end if 
!    end do

    do gp=1, ng-1
        if (xi<cdfsct(gp, g, tcomp)) exit
    enddo
    
    sample_scattering_group = gp
end function

subroutine sample_Nneutron(k, ke, nn, nnw, nwt)
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    real(8), intent(in) :: k, ke
    integer, intent(out) :: nn, nnw
    real(8) ::x_T, xi, nuf
    real(8), intent(out) :: nwt
    
    x_T = Xtot(tcomp, g)
    nuf = Xnuf(tcomp, g)
#ifdef _OPENMP
    if(unifsrc) then
        nwt=pwadjr(lc, mc)
        nn = pwadj(lc, mc)/k * nuf/x_T*w + GetRNS(lseed)
        nnw=0
!        nnw = 1._8/ke * nuf/x_T*w*fiswtmap(lc,mc) + GetRNS(lseed)
    else
        nwt=1.0_8
        nn = 1._8/k * nuf/x_T*w + GetRNS(lseed)
        nnw=0
        testnnt=testnnt+nn
!        nnw = 1._8/ke * nuf/x_T*w + GetRNS(lseed)
    endif
    tlk = tlk + nuf/x_T*w
    
#else
    nn = 1._8/k * nuf/x_T*w + GetRN()
    nnw = 1._8/ke * nuf/x_T*w + GetRN()
    gk = gk + nuf/x_T*w
#endif    
!    nfsrcth(lc,mc,tid)=nfsrcth(lc,mc,tid)+nn
end subroutine

subroutine addntrn(iq)
    include "pspec.FI"
    include "mcbase.FI" 
    integer, intent(in) :: iq
    real(8) :: a,ms,mu
    integer :: i, j, n
#ifdef _OPENMP
    n = tlnq(iq) +1
    nfsrcth(lc,mc,tid)=nfsrcth(lc,mc,tid)+1
    nfsrcfth(l,m,tid)=nfsrcfth(l,m,tid)+1
!    n = tlnq(iq) +1
    tlntrn((n),iq)%xy = xy
    tlntrn((n),iq)%w = nwt
    tlntrn((n),iq)%lm = lm
    tlnq(iq)=n
!    tlntrn((n),iq)%g = nn_group()
!    mu = 1-2*GetRNS(lseed)
!    ms = sqrt(1._8-mu**2)
!    a = 2*pi*GetRNS(lseed)
!    tlntrn((n),iq)%mu  = mu
!    tlntrn((n),iq)%dir(1) = ms*cos(a)
!    tlntrn((n),iq)%dir(2) = ms*sin(a)
!    tlnq(iq) = n  
!!    tlntrn((n),iq)%seed = stride(lseed)
!    call long2short(stride(lseed),tlntrn((n),iq)%seed)
!!    tlnq(iq) = n  
#else    
    n = nq(iq) +1 
    ntrn((n),iq)%xy = xy 
    ntrn((n),iq)%w = 1. 
    ntrn((n),iq)%lm = lm 
!    ntrn((n),iq)%g = nn_group()
!    mu  = 1-2*GetRN()
!    ms = sqrt(1._8-mu**2)
!    a = 2*pi*GetRN()
!    ntrn((n),iq)%mu  = mu
!    ntrn((n),iq)%dir(1) = ms*cos(a)
!    ntrn((n),iq)%dir(2) = ms*sin(a)
!    call long2short(stride(lseed),tlntrn((n),iq)%seed)
!    nq(iq) = n  
#endif     
end subroutine

integer function nn_group()
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    integer :: g
    real(8) :: sum, xi
#ifdef _OPENMP    
    xi  = GetRNS(lseed)
#else
    xi  = GetRN()
#endif
    sum = 0.
    do g=1, ng-1
        sum = sum + Xchi(tcomp,g)
        if (xi<sum) then
            exit
        end if 
    end do
    
    nn_group = g
end function    

end subroutine  
!End of Migration
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

! Event -> surface crossing
subroutine TFMC_SURF(surf, dir, tcomp, g, lm, x, y, w)
    implicit none
    include "FMC.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(in) :: surf, tcomp, g, lm(2)
    real(8), intent(in) :: dir(2), x, y, w
    integer :: i, d, n, lmo
    real(8) :: phs, od, on, th
    integer :: lmn(2)
    real(8) :: deno, h(2)
    integer :: signx(2,4)
    data signx /+1,-1, -1,+1, -1,+1, +1,-1/
    
    h(1) = hfmx
    h(2) = hfmy
    
    d = LTS2D(surf)
    n = LTS2N(surf)
    od = dir(d)
    on = dir(n)
    deno = 1._8/xtot(tcomp,g)/vfm
    phs = 1._8/abs(od)*deno

    do i=1, 2
        lmn = lm + LTFST(:,i,surf)
        if(d .eq. 1) then ! x-dir
            th = abs(hacy(lmn(2)-1)-y)
        else
            th = abs(hacx(lmn(1)-1)-x)
        endif
        tphi2(d,lmn(1),lmn(2)) = tphi2(d,lmn(1),lmn(2))+th*phs*od*od
        tj(d,lmn(1),lmn(2)) = tj(d,lmn(1),lmn(2))+ xtot(tcomp,g)*h(d)*phs*abs(od)
        txsurf(d,lmn(1),lmn(2)) = txsurf(d,lmn(1),lmn(2)) + h(d)*phs*od*on*signx(i,surf)
    enddo
end subroutine

! Event -> in mesh
subroutine TFMC_INMESH(a,b,x0,x1,trkl, dir, lm, g, w, tcomp)
    implicit none
    include "FMC.FI"
    include "pspec.FI"
    include "xsec.FI"
    real(8),intent(in) :: a, b, x0, x1, trkl, dir(2), w
    integer,intent(in) :: lm(2), tcomp, g
    integer :: surf
    integer :: i, j, d, lmo(2), lmn(2), l, m
    real(8) :: phs, oe, oo, th
    real(8) :: ap, bt, gm, xp, yp, ps
    real(8) :: mid(2), xx(4),xxx(4), xm, ym, ox
    integer :: cnt, lmFMC(2)
    
    l = lm(1)
    m = lm(2)
    ox=dir(1)*dir(2)
    do i=1, 4   
        lmo = lm + LTITER(:,i)
        lmn = lm + LTMI(:,i)
!        tlfmc(FXA,i,l,m) = tlfmc(FXA,i,l,m) + trkl*dir(1)*dir(2)
        xp = hacx(lmo(1))
        yp = hacy(lmo(2)) 
        ap = a/3
        bt = (b-yp-a*xp)/2
        gm = xp*yp-b*xp
        ps = abs((ap*(x1**3-x0**3) + bt*(x1**2-x0**2) + gm*(x1-x0))/vfm/dir(1)*w) !*sqrt(1-mu**2))
        tlfmc(FFA,i,l,m)  =tlfmc(FFA,i,l,m)  + ps
        tAbs(lmn(1),lmn(2))=tAbs(lmn(1),lmn(2)) + ps*xrmv(tcomp,g)
        tNuf(lmn(1),lmn(2))=tNuf(lmn(1),lmn(2)) + ps*xnuf(tcomp,g)
        txnavg(lmn(1),lmn(2))=txnavg(lmn(1),lmn(2)) + ps*ox*LTSIGNX(i)
    end do
    
    mid(1) = hacx(l)-hfmx/2
    mid(2) = hacy(m)-hfmy/2
    xx(1) = mid(1)
    xx(2) = (mid(2)-b)/a
    cnt = 0
    xxx = 0
    if (x0<x1) then
        xxx(inc(cnt))=x0
        if (xx(1)>x0 .and. xx(1)<x1) xxx(inc(cnt)) = xx(1)
        if (xx(2)>x0 .and. xx(2)<x1) then
            if (cnt .eq. 2 .and. xx(1)>xx(2)) then
                xxx(inc(cnt))=xx(1);xxx(cnt-1)=xx(2)
            else
                xxx(inc(cnt))=xx(2)
            endif
        endif
        xxx(inc(cnt))=x1
    else
        xxx(inc(cnt))=x1
        if (xx(1)<x0 .and. xx(1)>x1) xxx(inc(cnt)) = xx(1)
        if (xx(2)<x0 .and. xx(2)>x1) then
            if (cnt .eq. 2 .and. xx(1)>xx(2)) then
                xxx(inc(cnt))=xx(1);xxx(cnt-1)=xx(2)
            else
                xxx(inc(cnt))=xx(2)
            endif
        endif
        xxx(inc(cnt))=x0
    endif 
    do j=1, cnt-1
        xm = (xxx(j)+xxx(j+1))/2
        ym = a*xm+b
        lmFMC = lm + LTMI(:,mi(xm,ym))
        phi0(lmFMC(1),lmFMC(2))=phi0(lmFMC(1),lmFMC(2))+(xxx(j+1)-xxx(j))/abs(dir(1))*w*deno(lmFMC)
    end do
contains
function mi(xm, ym)
    integer :: mi
    real(8) :: xm, ym
    mi = 1
    if (xm>mid(1)) mi = mi+1
    if (ym>mid(2)) mi = mi+2
end function
function deno(lmFMC)
    real(8) :: deno
    integer :: lmFMC(2)
    deno = 1._8/vfm
    if (lmFMC(1).eq.1 .or. lmFMC(1).eq.nfmx+1) deno = deno*2
    if (lmFMC(2).eq.1 .or. lmFMC(2).eq.nfmy+1) deno = deno*2
end function
end subroutine



integer function nn_group_1(tc, xi)
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    integer :: g,tc
    real(8) :: sum, xi
    sum = 0.
    do g=1, ng-1
        sum = sum + Xchi(tc,g)
        if (xi<sum) then
            exit
        end if 
    end do
    
    nn_group_1 = g
end function    


subroutine switch_cycle(bat)
    integer :: bat(3)
    integer :: temp
    temp = bat(1)
    bat(1) = bat(2)
    bat(2) = temp
end subroutine
subroutine switch_wielandt(bat)
    integer :: bat(3)
    integer :: temp
    temp = bat(1)
    bat(1) = bat(3)
    bat(3) = temp    
end subroutine
function shiftk(iter, delk) result (ret)
    real(8) :: SHIFT_GUIDe(18)
    data SHIFT_GUIDE /1.0_8, 1.0_8, 0.75_8, 0.50_8, 0.35_8, 0.25_8, &
                    0.20_8, 0.15_8, 0.125_8, 0.100_8, 0.075_8, 0.050_8, &
                    0.035_8, 0.025_8, 0.020_8, 0.015_8, 0.0125_8, 0.0100_8/
    integer, intent(in) :: iter
    real(8), intent(in) :: delk
    real(8) :: ret
    ret = delk
    if (iter .le. 18) ret = SHIFT_GUIDE(iter) 
end function

function ShannonEntropy(ndimx, ndimy, psif)
    implicit none
!    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    real(8) :: ShannonEntropy, pj, fsum
    real(8), pointer, dimension(:,:) :: psif
    integer :: i, j, g, ndimx, ndimy
	ShannonEntropy = 0.
	fsum = sum(psif)
	do j=1, ndimy
	    do i=1, ndimx
    	    pj =psif(i,j)/fsum
    	    if (pj>0.) ShannonEntropy = ShannonEntropy -pj*log(pj)
    	end do
	end do
	ShannonEntropy = ShannonEntropy/ log(2._8)
end function


subroutine initp_from_dist(nht, ntrn)
#ifdef MPIF
	use MPIComm_mod
#endif
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    integer, intent(inout) :: nht
    type(particle), dimension(:), intent(inout) :: ntrn
    integer :: n, i, j, l, nq, en, tpin, tcomp
    real(8) :: wt, a, ms, mu
    integer(8) :: bseed, seed
    integer :: nwloc, left, npstride 
    real(8) :: xy(2), xybase(2)
    type (strcpin), pointer :: ppin 

    bseed=gseed
    nwloc = nht
#ifdef MPIF ! if MPI is used, seeds will be skipped. 
	nwloc = nht/nprocs
	left = nht-nprocs*nwloc
	nqmpi=nwloc
	nqmpi(0:left-1)=nwloc+1
	if (left >myrank) then 
	    nwloc=nwloc+1
	    npstride=nwloc*myrank
	    bseed=strideN(bseed, npstride)
	else
	    npstride=nwloc*myrank+left
	    bseed=strideN(bseed, npstride)
	endif
#endif
    nq = 0
    do l=1, nwloc
        seed=bseed
        nq = nq+1
        tpin = 0
        do while(tpin .eq. 0)
            i = nfmx*getRNS(seed)+1
            j = nfmy*getRNS(seed)+1
            tpin = cplst(i,j)
            if (tpin>0) then
                tcomp = typepin(tpin)%rc(typepin(tpin)%nr)
                if (Xnuf(tcomp,1) .eq. 0) tpin=0
            endif
        enddo
        ppin=>typepin(tpin)
        if (ppin%nr .eq. 0) then
            xy(1)=hfmx*(GetRNS(seed)-0.5_8)
            xy(2)=hfmy*(GetRNS(seed)-0.5_8)
        else
            call samplerxy(ppin%rr(ppin%nr), xy(1), xy(2), seed)
        endif
        xybase(1)=hacx(i-1)+hfmx/2
        xybase(2)=hacy(j-1)+hfmy/2
        ntrn(nq)%xy=xybase+xy
        ntrn(nq)%w = 1._8
        ntrn(nq)%lm(1) = i
        ntrn(nq)%lm(2) = j
        bseed=stride(bseed)
    end do
    nht = nq
end subroutine  

subroutine FS_FDM(nqn, nqtot, nwsum, ntrn) !======================================================================
#ifdef MPIF
    use MPIComm_mod
#endif
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(inout) :: nqn, nqtot
    integer, intent(in) :: nwsum
    type(particle), dimension(:), intent(inout) :: ntrn
    integer :: i,j, ncount, left
    integer :: noffset ! used when MPI is not applied
    integer(8) :: seed ! local seed for generating MC FS
    real(8) :: expfs ! expected # of FS
!    integer (8) :: bseed, seed, sdfs
    integer :: xb, yb
    integer :: nfsr(nfmx) !log for initial fission source distribution
#ifdef MPIF
!    integer :: nwloc, nqloc(0:nprocs-1), state
    integer :: itag, ni
    integer, dimension(0:nprocs) :: taginc
 
    do i=0, nprocs-1
        taginc(i)=i+1
    enddo
    taginc(i)=0
    itag=0
    nqmpi=0

    call bcast(initFSD, nfmx, nfmy, mpi_mc_comm)
#endif

    seed=gseed
    nqn = nwsum
    noffset=0
    open (71, file="nfsr.txt", status='unknown')
    do j=1, nfmy
        nfsr=0
        do i=1, nfmx
            expfs=nqn*initFSD(i,j)
            ncount=int4(expfs)
            if((expfs-ncount)>GetRNS(seed)) ncount=ncount+1
#ifdef MPIF
            do ni=1, ncount
                call addN2locQueue
            enddo
#else
            if(ncount>0) &
                call resampling(i,j,ntrn,noffset,ncount, seed)
            noffset=noffset+ncount
#endif
            nfsr(i)=cplst(i,j)
        end do
        write(71, '(2000(i10,1x))'), nfsr
    end do
    close(71)
#ifdef MPIF
    nqn=nqmpi(myrank)
    nqtot=sum(nqmpi)
    call init_weight(nqn, nwsum, ntrn)
    contains
    subroutine addN2locQueue
        if (myrank .eq. itag) then
            call resampling(i,j,ntrn,nqmpi(myrank),1,seed)
        else
            seed=stride(seed)
        endif
        nqmpi(itag)=nqmpi(itag)+1
        itag=taginc(itag)
    endsubroutine
#else
    nqn=noffset
    nqtot=nqn
    call init_weight(nqn, nwsum, ntrn)
#endif

end subroutine !======================================================================

subroutine FS_CMFD(nqn, nqtot, nwsum, ntrn) !======================================================================
#ifdef MPIF
	use MPIComm_mod
#endif
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(inout) :: nqn, nqtot
    integer, intent(in) :: nwsum
    type(particle), dimension(:), intent(inout) :: ntrn
    integer :: i,j,noffset, ncount, left
    real(8) :: expfs
    integer :: xb, yb
    integer(8) :: seed
	nqn = nwsum
    noffset=0
!    bseed=getSEED()
    seed=gseed

#ifdef MPIF
	call bcast(adjp, nfmx, nfmy, mpi_mc_comm)
	nqn = nwsum/nprocs
	left = nwsum-nprocs*nqn
	xb=myrank+1
	yb=nprocs-xb
	if (left >myrank) then 
	    nqn=nqn+1
	    bseed=strideN(bseed, nqn*myrank)
	else
	    bseed=strideN(bseed, nqn*myrank+left)
	endif

    do j=yb, nfmy, nprocs
        do i=xb, nfmx,nprocs
            expfs=nqn*adjp(i,j)
            ncount=int4(expfs)
            if((expfs-ncount)>GetRN()) ncount=ncount+1
            if(ncount>0) call resampling(i,j,ntrn,noffset,ncount, seed)
            noffset=noffset+ncount
        end do
    end do
    call nbalance(ntrn, nqn, nqtot)
#else
    do j=1, nfmy
        do i=1, nfmx
            expfs=nqn*adjp(i,j)
            ncount=int4(expfs)
            if((expfs-ncount)>GetRN()) ncount=ncount+1
            if(ncount>0) call resampling(i,j,ntrn,noffset,ncount, seed)
            noffset=noffset+ncount
        end do
    end do
    nqn=noffset
#endif

    
    ! population control
    call init_weight(nqn, nwsum, ntrn)
end subroutine !======================================================================

subroutine resampling(mi,mj,ntrn, noffset,ncount, bseed) !======================================================================
    include "pspec.FI"
    integer,intent(in) :: mi,mj,noffset, ncount
    real(8) :: xy(2), mu, ms, a
    integer :: i, tpin
    type(particle), dimension(:), intent(inout) :: ntrn
    type (strcpin), pointer :: ppin
    real(8) :: xybase(2)
    integer(8), intent(inout) :: bseed
    integer(8) :: seed
    xybase(1)=hacx(mi)-hfmx/2
    xybase(2)=hacy(mj)-hfmy/2
    ppin=>typepin(cplst(mi,mj))
    do i=noffset+1, noffset+ncount
        ! assuming that the inner most cell is a fuel pellet
        seed= bseed
        if (ppin%nr .eq. 0) then
            xy(1)=hfmx*(GetRNS(seed)-0.5_8)
            xy(2)=hfmy*(GetRNS(seed)-0.5_8)
        else
            call samplerxy(ppin%rr(ppin%nr), xy(1), xy(2), seed)
        endif
        mu = 1-2*GetRNS(seed)
        ms = sqrt(1._8-mu**2)
        a = 2*pi*GetRNS(seed)

        ntrn(i)%xy      = xybase+xy                                 
        ntrn(i)%w       = 1._8                                         
        ntrn(i)%lm(1)   = mi
        ntrn(i)%lm(2)   = mj
        bseed=stride(bseed)
    enddo
endsubroutine !======================================================================

subroutine samplerxy(r,x,y, seed)
    real(8), intent(in) :: r
    real(8), intent(out) :: x,y
    real(8) :: rsq
    integer(8), intent(inout) :: seed
    rsq=r**2
    ! 75.8% of acceptance
    do while(.true.)
        x=r*(1-2*GetRNS(seed))
        y=r*(1-2*GetRNS(seed))
        if (x**2+y**2<rsq) exit
    enddo
endsubroutine

subroutine LEAK2NETJ(ntallywt)
    include "pspec.FI"
    include "xsec.FI"
    include "arrays.FI"
    include "accm.FI"
    include "FMC.FI"
    integer, intent(in) :: ntallywt
    integer :: i, j, k, l, oi, oj, g
    real(8) :: hfm(4), nrm(4)
    real(8), dimension(4,nfmx, nfmy) :: jtemp
    real(8), allocatable, dimension(:,:,:,:) :: jcompmp
!! CHECK THE CONDENSED CURRENT =========================
!!    PRINT*, SUM(JF), SUM(HMJC)
!!    DO G=1, 23;DO J=1, NFMY;DO I=1, NFMX;DO K=1, 4
!!        HMJC(K,1,I,J)=HMJC(K,1,I,J)-jf(K,I,J,G)
!!    ENDDO;ENDDO;ENDDO;ENDDO
!!    PRINT*, SUM(HMJC(:,1,:,:))
!    PRINT*, SUM(JF), SUM(HMJC)
!    I=17;K=DE;
!    PRINT*, HMJC(K,1,1,1), SUM(jf(K,i,1:17,1:23))
!    I=18;K=DW;
!    PRINT*, HMJC(K,1,2,1), SUM(jf(K,i,1:17,1:23))
!    I=34;K=DE;
!    PRINT*, HMJC(K,1,2,1), SUM(jf(K,i,1:17,1:23))
!    I=34;K=DE;
!    PRINT*, HMJC(K,1,2,2), SUM(jf(K,i,18:34,1:23))
!
!    I=17;K=DE;
!    PRINT*, HMJC(K,2,1,1), SUM(jf(K,i,1:17,24:47))
!    I=18;K=DW;
!    PRINT*, HMJC(K,2,2,1), SUM(jf(K,i,1:17,24:47))
!    I=34;K=DE;
!    PRINT*, HMJC(K,2,2,1), SUM(jf(K,i,1:17,24:47))
!    I=34;K=DE;
!    PRINT*, HMJC(K,2,2,2), SUM(jf(K,i,18:34,24:47))
!!======================================================
    nrm(1:2) = 1._8/ntallywt/hcmx
    nrm(3:4) = 1._8/ntallywt/hcmy

!    do j=1, ncmy;do i=1, ncmx;do k=1, 4
!        oi = i + LTNEIGH(1,k)
!        oj = j + LTNEIGH(2,k)
!        hmjc(k,1,i,j) = (jcomp(k,1,i,j)-jcomp(LTOPST(k),1,oi,oj))*nrm(k)
!        hmjc(k,2,i,j) = (jcomp(k,2,i,j)-jcomp(LTOPST(k),2,oi,oj))*nrm(k)
!    enddo;enddo;enddo
    
    
    hfm(1:2) = hfmx
    hfm(3:4) = hfmy
    do g=1, ng
        do j=1, nfmy
            do i=1, nfmx
                do k=1, 4
                    oi = i + LTNEIGH(1,k)
                    oj = j + LTNEIGH(2,k)
                    jtemp(k,i,j) = (jf(k,i,j,g)-jf(LTOPST(k),oi,oj,g))/hfm(k)
                    
!                    call neighf2(i,j,k,oi,oj)
!                    if (bndryf(oi, oj).eq. 0) then
!                        jtemp(k,i,j) = (jf(k,i,j,g)-jf(LTOPST(k),oi,oj,g))/hfm(k)
!                    else
!                        jtemp(k,i,j) = jf(k,i,j,g)/hfm(k)
!                    end if
#ifdef _FMC
                    if (bndryf(oi, oj).eq. 0) then
                        do l=1, 8
                            ffmc(l,k,i,j) = (tlfmc(l,k,i,j)+tlfmc(l,LTOPST(k),oi,oj))/ntallywt/hfm(k)
                        enddo
                    else
                        do l=1, 8
                            ffmc(l,k,i,j) = tlfmc(l,k,i,j)/ntallywt/hfm(k)
                        enddo
                    endif
                    do l=9,12
                        ffmc(l,k,i,j) = tlfmc(l,k,i,j)/ntallywt/hfm(k)
                    enddo
#endif            
                end do
            end do
        end do
        do j=1, nfmy
            do i=1, nfmx
                jf(:,i,j,g)=jtemp(:,i,j)
                jf(:,i,j,g) = jf(:,i,j,g) / ntallywt
                jfAC(:,i,j,g) = jfAC(:,i,j,g) + jf(:,i,j,g)
#ifdef _FMC
                tlfmc(:,:,i,j)=0._8
#endif
            end do
        end do
    end do
end subroutine

subroutine evaluate_flux(nht)
    implicit none
    include "arrays.FI"
    include "accm.FI"
    include "pspec.FI"
    include "xsec.FI"
    include "FMC.FI"
    include "statistics.FI"
    integer, intent(in) :: nht
    integer :: g,i,j,l
    real(8) :: nrm, nrmc, rn
    real(8) :: fr(4)
    real(8) :: dif, dfsqsum(2), nrmf(2)
    integer :: hfx
    integer :: tcnt(4)
    integer :: tpin, nr
    real(8) :: vr, pinflux, rvpin, rvr
    autocor=0
    nrm = 1._8/nht
    psif = 0.
    ppf = 0.
    rvpin=1._8/vpin
    
    nrmc=1._8/nht/vcm
    do j=1, ncmy;do i=1, ncmx
        hmphi(1,i,j)=hmphit(1,i,j)*nrmc
        hmphi(2,i,j)=hmphit(2,i,j)*nrmc
    enddo;enddo
    
    do g = 1, ng
        do j=1, nfmy
            do i=1, nfmx
                tpin=cplst(i,j)
                if (tpin .eq. 0) cycle
                nr=typepin(tpin)%nr
                pinflux=0

                do l=0, nr
                    vr=typepin(tpin)%rv(l)
                    rvr=typepin(tpin)%rvr(l)
                    phif(l,i,j,g) = trk_est(l,i,j,g)*nrm*rvr
!                    pinflux = pinflux+phif(l,i,j,g)*vr 
                    phivac(l,i,j,g)=phivac(l,i,j,g)+phif(l,i,j,g)
                    psif(i,j) = psif(i,j) + Xnuf(typepin(tpin)%rc(l),g)*phif(l,i,j,g)*vr;
!                    ppf(i,j)  = ppf(i,j)  + Xkappa(typepin(tpin)%rc(l),g)*phif(l,i,j,g)*vr;
                enddo
!                stphi(i,j,g)   = stphi(i,j,g)+pinflux*rvpin
!                stphisq(i,j,g) = stphisq(i,j,g)+pinflux**2
           end do
        end do
    end do

    do j=1, nfmy
        do i=1, nfmx
            psifAC(i,j) = psifAC(i,j) + psif(i,j)
       end do
    end do

    do j=1, nfmy
        do i=1, nfmx
            stpsi(i,j)   = stpsi(i,j)   + psif(i,j)
            stpsisq(i,j) = stpsisq(i,j) + psif(i,j)**2
!            stpp(i,j)    = stpp(i,j)    + ppf(i,j)
!            stppsq(i,j)  = stppsq(i,j)  + ppf(i,j)**2
        end do
    enddo
!    write(fsptr, '(2000(e11.5,1x))') (stphi(1,1,i), i=1, ng)
end subroutine

subroutine estimateK (nht, k)
    include "pspec.FI"
    integer, intent(in)  :: nht
    real(8), intent(out) :: k
    k = gk / nht
end subroutine

subroutine cmfd_ctrl(lact, batch, oncmfd, onfdb) 
    implicit none
    include "pspec.FI"
    logical, intent(in) :: lact
    integer, intent(in) :: batch
    logical, intent(out) :: oncmfd, onfdb
    integer, save :: skip=0
    integer, save :: ninc=0
    oncmfd = .false.
    onfdb =  .false. 
!    skip = skip +1
!    if (skip .ge. nskip ) then
!        oncmfd = .true.
!        onfdb =  .true. 
!    endif
    
    if (lact .eq. .false. ) then
        if (CMFDIACT ) then
            skip = skip +1
            if (skip .ge. nskip)then
                oncmfd = .true.
                onfdb = .true.
                ninc = 0
            elseif (CMFDACT .and. skip .eq. inactive) then
                oncmfd =  .true.
                onfdb = .true.
            endif
            if (onfdb .and. batch<=ndelayFDB) onfdb = .false.
            if (fdbinact .eq. .false.) onfdb=.false.
        endif
    else
        if (CMFDACT) oncmfd = .true.
        if (fdbact) then
            onfdb = .true.
            ninc = ninc+1
            if (ninc .eq. nprd) then
                ninc = 0
                oncmfd = .true.
            endif
            if (RSTDHAT .and. onfdb .and. batch <=ndelayFDB) onfdb = .false.
        endif
    endif
endsubroutine

subroutine feedback(nqn, nwsum, ntrn)
#ifdef MPIF
	use MPIComm_mod
#endif
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(inout) :: nqn
    integer, intent(in) :: nwsum
    type(particle), dimension(:), intent(inout) :: ntrn
    integer ::i, j, ic, jc, lm(2), g, fi, fj
    real(8) w, wsum, ratio, wtest, wsumtot
    real(8) :: nsrc(ncmx, ncmy), nsrctot(ncmx, ncmy), wtadj(ncmx, ncmy)
    real(8) :: pj, shn, ms, a
    integer :: lc, mc
    real(8) :: psisum, rpsisum, upper, lower
    nsrc = 0
    do i=1, nqn
        lm = ntrn(i)%lm
        ic = f2cx(lm(1))
        jc = f2cy(lm(2))
        nsrc(ic,jc) = nsrc(ic,jc)+ntrn(i)%w
    end do
#ifdef MPIF
	call bcast(adjp, ncmx, ncmy, mpi_mc_comm)
	call reduce(nsrc, nsrctot, ncmx, ncmy, mpi_mc_comm, .true.)
	do j=1, ncmy;do i=1, ncmx
		nsrc(i,j)=nsrctot(i,j)
	enddo;enddo
!	if (myrank .eq. 0) print *, myrank, sum(adjp), sum(nsrctot)
#endif
    wtadj = 1
!    shn =0
    do j=1, ncmy
        do i=1, ncmx
            if (nsrc(i,j)>0) then
                wtadj(i,j) = adjp(i,j) * nwsum / nsrc(i,j)

!               pj= 1.0_8*nsrc(i,j)/nqn
!               shn = shn - pj * log(pj)
            elseif (adjp(i,j)>0) then
#ifdef MPIF
				if (myrank .eq.  nprocs-1) cycle ! add particles for the last process
#endif

                ! Add neutron at the coarse mesh
                fi = rcx *(i-1) + rcx*GetRN()+1
                fj = rcy *(j-1) + rcy*GetRN()+1
                nqn = nqn + 1
                ntrn(nqn)%xy(1) = hacx(fi-1)+hfmx*GetRN()
                ntrn(nqn)%xy(2) = hacy(fj-1)+hfmy*GetRN()
                ntrn(nqn)%w = 1._8                                         
                ntrn(nqn)%lm(1) = fi
                ntrn(nqn)%lm(2) = fj
!                ntrn(nqn)%g = nn_group_1(typepin(cplst(fi,fj))%rc(0), GetRN())
!                ntrn(nqn)%mu = 1-2*GetRN()
!                ms = sqrt(1._8-ntrn(nqn)%mu**2)
!                a = 2*pi*GetRN()
!                ntrn(nqn)%dir(1) = ms*cos(a)
!                ntrn(nqn)%dir(2) = ms*sin(a)
                
                wtadj(i,j) = adjp(i,j) * nwsum 
            endif
        end do
    end do
!    shn = shn/ log(2._8)
!    autocor(1) = shn
    wsum =0
    do i=1, nqn
        lm = ntrn(i)%lm
        ic = f2cx(lm(1))
        jc = f2cy(lm(2))
        ntrn(i)%w = ntrn(i)%w * wtadj(ic,jc)
        wsum = wsum + ntrn(i)%w
    end do
#ifdef MPIF
	call reduce(wsum, wsumtot, mpi_mc_comm, .true.)
	wsum=wsumtot
!	print *, myrank, wsumtot
#endif
    if (abs(wsum-nwsum)>1e-5_8) then
        ratio = 1._8*nwsum/wsum
        do i=1, nqn
            ntrn(i)%w = ntrn(i)%w * ratio
        end do
    endif

    ! check convergence of modifying factor
#ifndef MPIF
    do j=1, ncmy;do i=1, ncmx
        psicm(i,j,1)=0
    enddo;enddo
    do j=1, nfmy
        mc=f2cy(j)
        do i=1, nfmx
            lc=f2cx(i)
            psicm(lc, mc,1)=psicm(lc, mc,1)+psif(j,i)
            psisum =psisum+psif(j,i)
        enddo
    enddo
    rpsisum=1._8/psisum
    do j=1, ncmy;do i=1, ncmx
        psicm(i,j,1)=psicm(i,j,1)*rpsisum
    enddo;enddo
    
    ! comparisions
    upper=0;lower=0
    do j=1, ncmy;do i=1, ncmx
        if( psicm(i,j,1)>0 ) upper=upper+abs(adjp(i,j)/psicm(i,j,1))*adjp(i,j)
    enddo;enddo
#endif
    
endsubroutine

subroutine init_weight(nqn, nwsum, ntrn)
#ifdef MPIF
	use MPIComm_mod
#endif
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer, intent(in) :: nqn, nwsum
    type(particle), dimension(:), intent(inout) :: ntrn
    integer ::i, l, g
    real(8) :: w, wsum, ratio, wstot
    
    !w = 1.0_8*nwsum/nqn
    wsum = 0
    do i=1, nqn
        wsum=wsum+ntrn(i)%w
    end do
#ifdef MPIF
	call reduce(wsum, wstot, mpi_mc_comm, .true.)
	wsum=wstot
#endif
	w = 1.0_8*nwsum/wsum
    do i=1, nqn
        ntrn(i)%w = ntrn(i)%w*w
    end do
end subroutine

subroutine adjustWeight(wtsum, nHistory, ntrn, nptc, wtm)
    use define
    implicit none
    include "arrays.FI"
    include "pspec.FI"
    include "xsec.FI"
    real(8), intent(inout) :: wtsum
    integer, intent(inout) :: nHistory
    integer, intent(in) :: nptc
    real(8), intent(in) :: wtm
    type(particle), dimension(:), intent(inout) :: ntrn
    integer ::i
    
    wtsum = wtsum*wtm
    nHistory = nHistory*wtm
    
    do i=1, nptc
        ntrn(i)%w = ntrn(i)%w * wtm
    end do
end subroutine



subroutine TFMC_TRK
    include "pspec.FI"
end subroutine

subroutine source_balance(n,bank)
    integer :: n
    integer, dimension(:) :: bank
    real(8) :: mean
    integer :: sum, expt, over
    integer,dimension(n) :: bank_n, bank_d
    integer :: i
    do i=1, n
        sum = sum+bank(i) 
    end do
    mean = 1._8*sum / n
    expt = int4(mean)
    over = sum - expt*n
    bank_n = expt
    do i=1, over
        bank_n(i) = bank_n(i)+1
    end do
    do i=1, n
        bank_d(i) = bank(i) - bank_n(i)
    end do
    
    do i=2, n
        if (bank_d(i) > 0) then
            call move(bank, i, 1, bank_d(i))
        endif
    end do
    
    do i=2, n
        if (bank_d(i) < 0) then
            call move(bank, 1, i, -bank_d(i))
        endif
    enddo
end subroutine

subroutine move(bank, src, dest, cnt)
    integer, dimension(:) :: bank
    integer :: src, dest, cnt
    bank(src) = bank(src) - cnt
    bank(dest) = bank(dest) + cnt
end subroutine

subroutine write_result
    include "arrays.FI"
    include "xsec.FI"
    include "pspec.FI"
    include "statistics.FI"
    integer,parameter :: favg=901, fstdev=902
    integer :: g, j, i, l
    real(8) :: rnr
    character(256) :: sfile1, sfile2
    real(8), dimension(nfmx) :: stdevbuf
    integer :: nfpin
    real(8) :: fsum, fnrm
    rnr = 1._8

    open(1111,file ="p2b.txt",status='unknown')
    do j=nfmy, 1, -1
        write(1111, '(2000(e11.5,1x))') (pin2bf(i,j), i=1, nfmx)
    enddo

    open(favg,file ="psi.txt",status='unknown')
    open(fstdev,file ="stdevpsi.txt",status='unknown')
    nfpin = 0
    fsum = 0._8
    do j=nfmy, 1,-1
        ! for error of flux
        do i=1, nfmx
            if (stpsi(i,j)>0._8) then
                stdevbuf(i) = rnr*sqrt((stpsisq(i,j)-stpsi(i,j)**2/active)/active/(active-1))
                nfpin = nfpin+1 
                stpsi(i,j) = stpsi(i,j)*rnr/active
                stdevbuf(i) = stdevbuf(i)/stpsi(i,j)
            else
                stdevbuf(i) = 0._8
            endif
            
            fsum = fsum + stpsi(i,j)
        enddo  
        ! write
        write(fstdev, '(2000(e11.5,1x))') (stdevbuf(i), i=1, nfmx)
    enddo
    
    fnrm = nfpin/fsum
    do j=nfmy, 1,-1
        write(favg, '(2000(e11.5,1x))') (stpsi(i,j)*fnrm, i=1, nfmx)
    enddo
    close(favg)
    close(fstdev)
    
    
#ifdef ASMPRL  
    open(favg,file ="nfs.txt",status='unknown')
    do j=ncmy, 1,-1
        write(favg, '(2000(i17,1x))') nfs(:,j)
    enddo
    close(favg)

    open(favg,file ="nin.txt",status='unknown')
    do j=ncmy, 1,-1
        write(favg, '(2000(i17,1x))') nin(:,j)
    enddo
    close(favg)

    open(favg,file ="nout.txt",status='unknown')
    do j=ncmy, 1,-1
        write(favg, '(2000(i17,1x))') nout(:,j)
    enddo
    close(favg)
    
    open(favg,file ="nabsr.txt",status='unknown')
    do j=ncmy, 1,-1
        write(favg, '(2000(i17,1x))') nabsr(:,j)
    enddo
    close(favg)
    
    !print *, sum(nfs), sum(nabsr), sum(nout), sum(nin)
    !print *, sum(nfs)-sum(nabsr), sum(nout)-sum(nin)
#endif
    
    
    open(favg,file ="psiz.txt",status='unknown')
    open(fstdev,file ="stdevpsiz.txt",status='unknown')
    nfpin = 0
    fsum = 0._8
    do j=nfmy, 1,-1
        ! for error of flux
        if (stpsiz(j)>0._8) then
                stdevbuf(1) = rnr*sqrt((stpsizsq(j)-stpsiz(j)**2/active)/active/(active-1))
                nfpin = nfpin+1 
                stpsiz(j) = stpsiz(j)*rnr/active
                stdevbuf(1) = stdevbuf(1)/stpsiz(j)
        else
                stdevbuf(1) = 0._8
        endif
            
        fsum = fsum + stpsiz(j)
        write(fstdev, '(2000(e11.5,1x))') stdevbuf(1)
    enddo
    fnrm = nfpin/fsum
    do j=nfmy, 1,-1
        write(favg, '(2000(e11.5,1x))') stpsiz(j)*fnrm
    enddo
    close(favg)
    close(fstdev)


    open(favg,file ="psia.txt",status='unknown')
    open(fstdev,file ="stdevpsia.txt",status='unknown')
    nfpin = 0
    fsum = 0._8
    do j=ncmy, 1,-1
        ! for error of flux
        l=(j-1)*ncmx
        do i=1, ncmx
            l=l+1
            if (stpsia(l,1)>0._8) then
                stdevbuf(i) = rnr*sqrt((stpsiasq(l,1)-stpsia(l,1)**2/active)/active/(active-1))
                nfpin = nfpin+1 
                stpsia(l,1) = stpsia(l,1)*rnr/active
                stdevbuf(i) = stdevbuf(i)/stpsia(l,1)
            else
                stdevbuf(i) = 0._8
            endif
            
            fsum = fsum + stpsia(l,1)
        enddo
        ! write
        write(fstdev, '(2000(e11.5,1x))') (stdevbuf(i), i=1, ncmx)
    enddo
    fnrm = nfpin/fsum
    do j=ncmy, 1,-1
        write(favg, '(2000(e11.5,1x))') (stpsia(l,1)*fnrm, l=(j-1)*ncmx+1, j*ncmx)
    enddo

    
    close(favg)
    close(fstdev)
    if (ireft .eq. 2) then
        open (2001, file="REF.bin", form='binary', status='unknown')
        write(2001) stpsi*fnrm
        close(2001)

        open (2001, file="REF_P2BF.bin", form='binary', status='unknown')
        write(2001) pin2bf
        close(2001)
        
    endif

    open(favg,file ="pinpower.txt",status='unknown')
    open(fstdev,file ="stdevpinpower.txt",status='unknown')
    do j=nfmy, 1,-1
        ! for error of flux
        do i=1, nfmx
            if (stpp(i,j)>0._8) then
                stdevbuf(i) = rnr*sqrt((stppsq(i,j)-stpp(i,j)**2/active)/active/(active-1))
            else
                stdevbuf(i) = 0._8
            endif
            stpp(i,j) = stpp(i,j)*rnr/active
        enddo
    ! write
    write(favg, '(2000(e11.5,1x))') (stpp(i,j), i=1, nfmx)
    write(fstdev, '(2000(e11.5,1x))') (stdevbuf(i), i=1, nfmx)
    enddo

    close(favg)
    close(fstdev)

!	if (unifsrc) then
        open(favg,file ="nfsrc.txt",status='unknown')
        do j=ncmy, 1,-1
            write(favg, '(2000(i8,1x))') (nfsrc(i,j), i=1, ncmx)
        enddo

        close(favg)
        open(favg,file ="nfsrcf.txt",status='unknown')
        do j=nfmy, 1,-1
            write(favg, '(2000(i8,1x))') (nfsrcf(i,j), i=1, nfmx)
        enddo

        close(favg)
!	endif

end subroutine
subroutine evaluate_corcoef
    use ifport
    use util_STT
    include "pspec.FI"
    include "arrays.FI"
    integer :: i,j,k
    integer :: st, ed
    integer, parameter :: n=10, maxbuf=2000
    real(8) :: mean, Cr, sum2, rk, e1, e2, s2
    real(8) :: covlag, acrc(n-1), stds, stdt
    integer, save :: ibuf=0
    real(8) :: buf(maxbuf)
    if (ibuf>maxbuf) return

    ibuf=ibuf+1 
    buf(ibuf)=psif(nfmy,1)
    
    st = ibuf-n+1
    ed = st + n - 1
    if (st<1) return
    sum2 = 0
    do i=1, n-2
        Cr=0
        stds = stdevp(buf(:),st,ed-i)
        stdt = stdevp(buf(:),st+i,ed)
        covlag = covarlag(buf(:), st, ed, i)
        acrc(i) = covlag/(stds*stdt)
        sum2  = sum2 + (n-i)* covlag
    end do
    s2 = var(buf(:),st, ed)
    s2 = s2/n
    autocor(2) = sqrt(abs(2._8/n/(n-1)*sum2 / s2))
end subroutine


#ifdef MPIF

subroutine nbalance(ntrn, nq, nqtot)
    use MPIComm_mod, only : mpi_mc_comm
	use define
    implicit none
    include "mpif.h"
    include "pspec.FI"
    include "mpi_arrays.FI"
    
    type(particle), dimension(:), intent(inout) :: ntrn
    integer, intent(inout) :: nq
    integer, intent(out) :: nqtot
    integer :: ierr, i, maxn, minn, nqavg, nqlft, nqsend, nqadd
    integer, dimension(0:nprocs-1) :: rcvcnt, idisp, nadd, nbeg, nqrcv, nqideal
    real(8), parameter :: tolr = 0.1_8
    logical, save :: first=.true.

    ! gather the # of source neutrons for each process
    call mpi_allgather(nq, 1, mpi_integer, &
                    nqrcv, 1, mpi_integer, &
                    0, mpi_mc_comm, ierr)
    nqtot = sum(nqrcv)
    nqmpi=nqrcv

    if (myrank .eq. 0) then
        if (first) then
            open(779, file='nbalance.txt', status='unknown')
            first=.false.
        endif
        write(779,'(100(i12, 1x))') nqmpi(:)
    endif


    nqavg = nqtot/nprocs
    nqlft = nqtot - nqavg*nprocs
    nqideal=nqavg
    do i=1, nqlft
        nqideal(nprocs-i)=nqideal(nprocs-1)+1
    enddo
    nqmpi(1:nprocs)=nqideal(:)
    nqsend=max(0,nq-nqideal(myrank))

    ! set the # of neutron to send or receive
    maxn=0
    minn=history
    idisp(0)=0
    nbeg=0
    do i=0, nprocs-1
        if (nqrcv(i)>maxn) maxn=nqrcv(i)
        if (nqrcv(i)<minn) minn=nqrcv(i)
        nadd(i)=max(0, nqideal(i)-nqrcv(i))
        rcvcnt(i)=max(0, nqrcv(i)-nqideal(i))
        if (i>0) then
            nbeg(i)=nbeg(i-1)+nadd(i-1)
            idisp(i) = idisp(i-1)+rcvcnt(i-1)
        endif
    enddo
    if (myrank .eq. 0) then
        print '(100(i12, 1x))', nqrcv
        print '(100(i12, 1x))', nqideal
        print '(100(i12, 1x))', nadd
        print '(100(i12, 1x))', rcvcnt
        print '(100(i12, 1x))', nbeg
    endif
    if (.true.) then
        call mpi_allgatherv(ntrn(nq-nqsend+1), nqsend, t_ntrn, &
                            nbuf, rcvcnt, idisp, t_ntrn, &
                            mpi_mc_comm, ierr)

        ! assigning new source neutrons to the queue
        do i=1, nadd(myrank)
            ntrn(nq+i).xy=nbuf(nbeg(myrank)+i).xy
            ntrn(nq+i).w=nbuf(nbeg(myrank)+i).w
            ntrn(nq+i).lm=nbuf(nbeg(myrank)+i).lm
            !ntrn(nq+i)=nbuf(nbeg(myrank)+i)
        enddo
        nq=nq-nqsend+nadd(myrank)
        nq=nqideal(myrank)
    endif
    
end subroutine
#endif

subroutine testnn(ntrn,nq, nqtot)
	use define
	implicit none
	include "pspec.FI"
	include "mpi_arrays.FI"

    type(particle), dimension(:), intent(inout) :: ntrn
	integer, intent(inout) :: nq
	integer, intent(out) :: nqtot
!	ntrn(1)%w=1.1
	print *, myrank,  "ntrn test", nq
	print *, ntrn(1)%w


	nqtot =0

end subroutine


end module
